In this notebook we reproduce the 1D simulations presented by Morse and Nielsen in their seminal paper of 1971. In this paper they explore the development of the Weibel instability, by initializing a single species with a temperature anisotropy.
Specifically, a uniform density electron species is initialized with $u_{thx} = 0.05 \,\rm{c}$ and $u_{thy} = 0.25 \,\rm{c}$. We use the same simulation parameters as the paper, but increase the resolution by a factor of 3 to have more detail on the phase-space plots, which gives 150 cells and a cell size $\Delta x = 0.1 \,c \omega_p^{-1}$, with $3 \times 10^4$ particles.
Also note that we
In [2]:
import em1ds as zpic
import numpy as np
import matplotlib.pyplot as plt
# Thermal velocity
uth = [0.05,0.25,0.0]
electrons = zpic.Species( "electrons", -1.0, 200, uth = uth )
sim = zpic.Simulation( 150, box = 15.0, dt = 0.08, species = electrons )
The plots in the original paper distiguish particles in 3 groups according to their initial $u_y$ value: the first group holds the 1/4 of particles that is closest to the minimum $u_y$ injected, the second group holds the 1/2 of particles that is closest to $u_y = 0$ and the third group holds the 1/4 of particles that is closest to the maximum $u_y$ injected. To reproduce this, we turn off particle sorting (so the particles will retain their initial indexes) and store the indexes of the particles belonging to each of the groups:
In [2]:
# Disables sorting
electrons.n_sort = 0
# Selects particle groups for visualization
uy = electrons.particles['uy']
lim = uth[1] * 0.67449
idxa = uy <= -lim
idxb = (uy > -lim ) & (uy < lim )
idxc = uy >= lim
For simplicity we define a routine to generate a side-by-side plot of the particle phasespace, with particles divided into 3 groups as described above, and the magnetic field along the z direction, $B_z$, like the plots in Fig. 3 of the paper:
In [3]:
# Routine for generating plots
def vis( sim ):
f, (ax1,ax2) = plt.subplots(ncols = 2, sharey=True)
f.set_size_inches(8,6)
# Phasespace plot
x = ( electrons.particles['ix'] + electrons.particles['x']) * electrons.dx
ux = electrons.particles['ux']
ax1.plot( ux[idxa] - 1, x[idxa], '.', ms=1,alpha=0.4, label = "$u_x^L - 1$")
ax1.plot( ux[idxb] , x[idxb], '.', ms=1,alpha=0.4, label = "$u_x^C$")
ax1.plot( ux[idxc] + 1, x[idxc], '.', ms=1,alpha=0.4, label = "$u_x^R + 1$")
ax1.set_xlabel("$u_x$ [$m_e c$]")
ax1.set_title("$u_x - x$ phasespace")
ax1.legend()
ax1.grid(True)
ax1.set_ylabel("$x$ [$c\,\omega_n^{-1}$]")
# Magnetic field plot
ax2.plot( sim.emf.Bz, np.linspace(0, sim.box, num = sim.nx), label = "$B_z$")
ax2.set_xlim(left=-0.15,right=0.15)
ax2.grid(True)
ax2.set_xlabel("$B_Z$ field []")
ax2.set_title("Magnetic Field")
ax2.legend()
f.suptitle("t = {:g}".format(sim.t) + " $\omega_p^{-1}$")
plt.show()
We now run the simulation until the specified times and visualize results for each one:
In [4]:
sim.run(30.0)
vis(sim)
In [5]:
sim.run(40.0)
vis(sim)
In [6]:
sim.run(50.0)
vis(sim)
In [7]:
sim.run(100.0)
vis(sim)
In [8]:
sim.run(200.0)
vis(sim)
In [9]:
sim.run(400.0)
vis(sim)
In [6]:
electrons = zpic.Species( "electrons", -1.0, 1000, uth = [0.05,0.25,0.0] )
sim = zpic.Simulation( 150, box = 15.0, dt = 0.08, species = electrons )
To store the field energy we run the simulation with a customized loop to store these values at every time-step:
In [7]:
import math
tmax = 400
niter = int(math.ceil(tmax / sim.dt))
EneB = np.zeros(niter)
EneE = np.zeros(niter)
norm = 0.5 * sim.emf.nx * sim.box / sim.nx
print("\nRunning simulation up to t = {:g} ...".format(tmax))
while sim.t < tmax:
print('n = {:d}, t = {:g}'.format(sim.n,sim.t), end = '\r')
EneB[sim.n] = np.sum(sim.emf.Bx**2+sim.emf.By**2+sim.emf.Bz**2) * norm
EneE[sim.n] = np.sum(sim.emf.Ex**2+sim.emf.Ey**2+sim.emf.Ez**2) * norm
sim.iter()
print("\nDone.")
In [ ]:
plt.plot(np.linspace(0, sim.t, num = niter),EneB, label = "MAGNETIC")
plt.plot(np.linspace(0, sim.t, num = niter),EneE, label = "ELECTRIC")
plt.yscale('log')
plt.ylim(ymin=0.01)
plt.grid(True)
plt.xlabel("$t$ [$1/\omega_n$]")
plt.ylabel("Field energy [$m_e c^2$]")
plt.title("Electro-Magnetic field energy")
plt.legend()
plt.show()
As shown in the paper, the energy lies almost exclusively in the magnetic field. The plot also clearly shows the linear stage of the instability as well as the saturation fase. We see mulitple local maxima in the magnetic field energy that are associated with the coalescence of the trapped particles. However, we cannot fully reproduce the multiple local maxima shown in the original paper, nor the continuous growth in energy following the saturation of the instability.
In [ ]: